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Abstract. The Bertaut approach associated with charge spreading so as to 
enhance the rate of convergence of Coulomb series in crystals is extended to the 
case of an arbitrary multiple spreading with a given initial spreading function. 
It is shown that the effect of spreading may in general be treated as a uniform 
transformation of space, providing that zero mean potential as a universal spatial 
property is sustained. As a result, electrostatic potentials driven by different 
orders of multiple spreading can be obtained from the same energy functional in 
a consistent manner. It is found that the effect of multiple spreading gives rise 
to more advanced forms described, for example, by simple exponential decrease, 
but the functional description based on a Gaussian spreading turns out to be 
invariant. In addition, the effects of a multiple charge spreading based on 
simple exponential and Gaussian spreading functions arc compared as typical 
of molecular calculations. 



PACS numbers: 02.30.Lt, 02.30.Uu, 61.50.Ah, 61.50.Lt 
1. Introduction 

The problem of lattice summation of Coulomb series over crystal structures is principal 
for describing solid state. Apart from a large number of traditional approaches 
to this subject [l] El E], many novel proposals for solving this problem still arise 
HISlElEKHliniliniinilll]- Nevertheless, the classical Ewald approach [13] remains 
one of the most effective and so widespread [TJl [151 UHl HZj- This is the reason 
that understanding the nature of this efficiency is of great importance [18]. In 
particular, its relation to the effect of screening Coulomb potentials was revealed by 
Nijboer and De Wette [19]. As a result, some generalizations associated with different 
types of screening have been proposed [TBI EQ] ■ Another fruitful explanation of the 
foregoing efficiency is based on the idea of charge spreading proposed first by Ewald 
in his original paper I13| and developed further by Bertaut '21'. In particular, this 
generalized approach turns out to be expedient in applications to molecular dynamics 
[TBI l22] . Here we will discuss this treatment in more amount of detail so as to coincide 
known variations inherent in its implementation. 

In the original paper of Bertaut [21] the pair-wise Coulomb interaction is 
discussed. As a result, the double charge spreading naturally arises in that approach. 
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In particular, this feature results in the fact that the square of the Fourier transform 
describing the spreading function takes place in the sum over the reciprocal lattice 
contributing to the Coulomb energy [5T1 [25]. On the other hand, it turns out 
that a single charge spreading is still sufficient if the electrostatic potential is first 
considered [HI [25l [26] . As a result, the Fourier transform of the spreading function, 
but not its square, arises in the sum over reciprocal lattice vectors contributing to the 
energy within such a treatment. This fact was the subject of discussion ^25l [27] [28]. 
Altogether, it was shown that both the treatments are quite correct and can eventually 
originate the description proposed by Bertaut. Nevertheless, the original treatment 
proposed by Bertaut appears to be symmetric with respect to both the set of charges 
generating potentials and the similar set of charges interacting with those potentials. It 
is instructive that the latter situation may be regarded as some uniform transformation 
of space [31 [22]. 

In the present paper we extend the concept of charge spreading and introduce 
the regular idea of a multiple spreading, bearing in mind that this effect can always 
be addressed to the charge distribution generating the potential field. On the other 
hand, such a standpoint is not obviously unique and therefore various other points of 
view are discussed. In particular, the idea of spreading as a uniform transformation of 
space is developed in a general form. In addition, in the case of a multiple spreading 
the universal character of the Ewald approach is recognized. The effect of a multiple 
charge spreading on the Coulomb interaction between a couple of objects neutral on 
average is also discussed. 

2. Preliminaries 

Let us consider a crystal described by a charge distribution pdr) subject to 
translational symmetry. On the other hand, we can also introduce a local charge 
distribution p(r) attributed to a unit cell and driven by the natural condition of 
electrical neutrality 



where the integration is carried out over a volume V occupied by p{r). Note that 
the unit-cell parallelepiped of volume v, constituted of three noncomplanar vectors 
of elementary translations of a Bravais lattice at hand, is assumed to be contained 
completely in V. The latter is essential if p(r) is spread beyond that parallelepiped 
volume |5] . Then the overall charge distribution in question can be written as 



where the summation over i is performed over sites specified by vectors Ri appropriate 
to the Bravais lattice of interest. It is evident that representation ([2]) for pc{r) is not 
unique due to an optional choice of p{r) [3]. Nevertheless, form ([2|) is subject to 
translational symmetry as well. This is the reason that Pc{r) defined by can be 
cast in terms of a series over reciprocal lattice vectors h: 




(1) 




(2) 




(3) 



h 
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where the prime on the summation sign impUes that the contribution of /i = is 
actually omitted, as follows from formula ([7]) derived later on. The structure factor 
F{h), by definition, is determined as [211 HE] 

F{h) = -[ pc{r)exp{-2TTihr)dr. (4) 

« JcoU 

Here the integration is carried out over the unit-cell parallelepiped mentioned above. 
Substituting ^ into ([4]), we obtain 

F{h)^-[ drS^ p{r ~ Ri)exp{~2TTih{r ~ R,)). (5) 

Bearing in mind that the integration over the unit-cell parallelepiped along with 
the summation over i is transformed into the integration over all space [29] reduced 
eventually to V intrinsic to p{r), relation ([5]) is converted into 

F(h)^-[ p(r)exp(~2nihr)dr. (6) 
V Jv 

It is important that 

Fih = 0) = 0, (7) 

in agreement with ([T|). 

The electrostatic potential exerted by charge distribution ([2]) at a reference point 
r is of the form 

J |ri-r| 

where the prime on the integral sign stands for missing a singular contribution of any 
point charge if it happens at r. If we now make use of Poisson's equation for the 
Green function 

Vl-^— = -A7rS{ri-r), (9) 
|ri - r\ 

where the differentiation is performed with respect to r and 6(r) is the Dirac delta 
function, we readily confirm from ^ that U{r) is subject to the conventional Poisson's 
equation of the form 

VlUir) - -iirp.ir). (10) 

On the other hand, on inserting ^ into ([5]), the result can be written as 

Uir)^[Y^'Fih)J'-^^^^^^^ (11) 

h 

where the square brackets decorated by the prime imply missing the same singular 
term mentioned in ([8]). It is easy to show that 

exp{2nihri) dri _ 1 



Substituting (fT^ into pT|) . we obtain 
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However, if we insert relation ^ directly into ([5]), the result is as follows: 

Uir)- fj£^\Y.Pi^^-RO- (14) 
If we go over to a new variable of integration r' = ri — Ri now, then 

wr^(r2^ (15) 

^ Jv m 

where 

R, = R, + r'-r, (16) 

the prime on the summation sign in p5|) means that the singular contribution 
associated with a point charge at r must be still excluded. The asterisk over the 
summation sign points to the fact that the summation over large Ri is not yet 
defined properly in formula (fT5|) so as to be consistent with the absence of the h ^ 
contribution in expression (|13p . 

This inconsistency is the essence of the conditional convergence of Coulomb series 
in crystals. It is especially pronounced in the particular case of point-charge lattices 
described by [23 

p(r-)-E*'^(^-^^)' (17) 

where the summation over j is carried out over point charges qj belonging to a unit 
cell, located at positions bj and governed by the condition 

in agreement with Substituting ([TT]) into we deduce 

F{h) = i ^ cxp{-2TTihbj). (19) 

3 

If relation p9| is now inserted into ([13]) and relation pT| is inserted into (|15p , where 
equation (|16p is taken into account, then we obtain 



h,j i j 



where 



= |fi, + 6, -r|. (21) 

Here the contribution of r = bj is supposed to be excluded in the first relation on 
the right-hand side of (|20p and the same contribution at i?i = is to be excluded 
in the second issue. Both of the expressions in (PO)) describe the same potential by 
definition, so that the singularity associated with the contribution of large Ri must 
be resolved as it is prescribed by exclusion of the h ~ term. However, even in this 
case the convergence of the sum over h is not fast for general r [3D] . The event of r 
at which a point charge exists is an exclusion, where a compensating term enhancing 
the rate of convergence arises [31] . 
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3. Multiple charge spreading 

In order to enhance the rate of convergence of the series mentioned above, we extend 
the treatment of Bertaut [21] and define a modified unit-cell charge distribution as 
follows: 

p(")(r) = Ja{\r - r,\)ai\n - r2\) ■ ■ ■ a(|r„_i - r„|) 

X p{rn)dri...drn, (22) 

where n identical functions (j{\r\) are introduced. These functions spread the actual 
charge at every point in a consecutive manner and are normalized by the condition 



J cr{\r\) dr ^ 4:77 J a{r)r'^ dr = 1, 



(23) 



where r = |r|. Like ([5]), the modified overall charge distribution in the crystal then 
takes the form 

P^"'(r) = E^^"'(^-^')- (24) 

i 

This value can in turn be cast in terms of the Fourier transforms associated with the 
reciprocal lattice vectors: 

pi'^Hr) = J2' F'-''\h)exp{2nihr), (25) 

h 

F^"\h)^- [ p(")(ri)exp(-27rih.ri)dri. (26) 

V Jccll 

Substituting ([M]) into (pS)) . we obtain 

F^'''\h) = - j dr exp{-2Trihr) J dri . . . dr„cr{\r - ri\) 

X (j(\ri - r2|) . . . a{\ Mrn). (27) 

If we here go over to new variables of integration 

r'n-^i ^ rn-i ~ Tn, r[^ri~r2, r' = r - ri, (28) 

keeping in mind that 

r = r' + r; + + +r„, (29) 
then it is easy to show that 

^'(")(/i) =F(/i)5"(/i). (30) 
Here the function S{h) is defined by the relation 

S{h) ^ J a{\r\)exp{~2TTihr)dr, (31) 

with the evident properties 

S{h) = S{-h), 5(0) = 1, (32) 



in agreement with formula p3|) . 

The modified electrostatic potential appropriate to (|24p is naturally equal to 

J \ri-r\ 
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Indeed, the substitution of ^ into ([55)1 yields Poisson's equation 

V2[7(")(r) =-47rp(")(ri) (34) 

associated with (|22p . Comparing ([33]) with ^ and taking relations (fT3|) and (l30|) into 
account, we readily derive 

C^^-Hr-) = -E^ ^^?.^2^^^ exp(27rzM, (35) 

^ h ' ' 

where any restriction associated with a point charge contribution is immaterial now 
due to the attenuation effect of S{h). This is a direct consequence of the fact that 
there are no point charges after transformation (j22p . 

On the other hand, if we insert definition (|22|) into (|33p and make use of relations 
(128)1 and dUl, then we obtain 



U 



'"^M=E / Pir')0'^"\\mdr', (36) 
^ Jv 

where definition ()16|1 is utilized, the asterisk over the summation sign points out that 
the problem of remote Hi still exists in the present representation, 
^H(i^l) ^ / ^(|n|)...a(|r„|)dn. .rfr-„^ 

To overcome the problem of remote the initial electrostatic potential of 
interest can be rewritten in the form 

[/(„) (r ) = ^(") (r) + [{/(r) - J/*") (r)] . (38) 

If we now utilize result ([55]) for the first term on the right-hand side of formula ([55)1 
and employ ([T5)) . ([TB)) and ([55)1 for the remainder, we get 

i l-Ril -''•=''3 

Here we introduce the following compact definition 

(40) 

for the difference characteristic of the case. As a result, the asterisk over the 
summation sign can be omitted, because the summation over i in (|39p is now carried 
out in a consistent manner resolving the conditional convergence of this sum at large 
Ri. On the other hand, the prime on the summation sign over i in (|39p stands for 
the omission of the singular contribution of a point charge, if it happens, provided 
that such a contribution would be described by the first term on the right-hand side 
of ([40)1 . Finally, the last term on the right-hand side of ()39p describes the elimination 
of the same contribution from the regular part specified by il^^^R) in ([40)1 . 

It is important that the convergence of the first term on the right-hand side of 
expression ([39| is expected to be rather fast due to the effect of S{h) and the same is 
right for the direct sum in the remainder, in accord with pTl, as will be discussed in 
more amount of detail in [55] . 

Note that in the case of n = 1 result ([5^ is tantamount to the Nijboer-De Wette 
approach |19| . bearing in mind that normalization (|23p is not principal here due to 
the fact that this representation is after all incorporated by means of identity ([38)l. 
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4. Charge spreading as a uniform transformation of space 

It is important that the spreading at hand can be treated as a uniform transformation 
of space [3]. Indeed, according to ([22|) . this transformation connecting an initial point 
r' with a final point r is of the form 



/(")(r,r') = J a{\r - r,\)a{\n - r^l) . . .a(|r(„_2) - nn^i)\) 

X (T(|r„_i-r'|)dri...dr„_i, (41) 
where the limiting cases of n = and n = 1 can be defined, respectively, as 

/(")(r,0=<5(r-r'), (42) 
/(i)(r,r') =CT(|r-r'|). (43) 
Note that in terms of (|¥T1) . definition takes the form 

(r) = j (r , r')p(r') dr' . (44) 

According to (j23|) . one can see that 

j f^'''\r,r')dr' (45) 

Moreover, it is evident from definition (|4ip that 

/(")(r,r') -/("H^',^)- (46) 
This symmetry implies that transformation (j4ip may be regarded either as a spreading 
of initial points containing charges or as a spreading of final points which may be free 
from charges. Furthermore, relation (|4ip can be represented as a following convolution: 



&\r,r') = j /(™)(r,ri)/("-™)(ri,r')dri, (47) 

where < to < n, with including the limiting cases specified by (|42| and (|43p in the 
integrand. 

Another important convolution arises from (|37p as connecting the initial point r' 
and the final point r in (fT6|) . Indeed, we can substitute expression (fT6|) in place of R 
in formula ([37|) and go over to the following new variables of integration : 

ri=r-r'i, r2=ri-r^, r„, ^ r'^^_^ - r'^, (48) 

where < m < n again. Integrating over r'^ and keeping relation (I4ip in mind, we get 

/2H(|H, +.'_.,). ^/'"-Hr-,n)/(— )(r.,rOdn^r-._ 
J |iii + r2-ri| 

Inserting (|47l) into (|44)) and ([50]) into ((36)) . one can readily show that 

(r) = y" (r , r')/5("-") (r') dr', (51) 

[7(")(r)== / /^™^(r,r')C/'"~"nr')'^^', (52) 



where < m < n and definitions (I44p and (j36p are, respectively, used in the integrands. 
Note that relations ([5T|) and ([5^ are of the same structure. Moreover, they turn out 
to be complementary to each other. The latter fact becomes evident if we consider 
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the bulk Coulomb energy per unit cell, which can be written down in a traditional 
fashion [551 [Ml 133] as 

£in)^\j^p{r)U^n){r)dr. (53) 

On substituting ([5^ into ([53)1 and taking equation ([5]) into account, relation ([5 
easily converted into 

\F{h)\^S''(h) 



h ' ' 



J^drdr'p{r)p{r') 
ql (54) 



^\ 

where the last term describes the correcting contribution of all point charges in the 
unit cell. Upon investigating the first term on the right-hand side, we may notice that 
the numerator of the summand can be represented in the form 

|F(/i)|2S'"(/i) ^ /'('")(/i)F("-'")*(h.), (55) 

where < to < n, in agreement with ([50)1 . In other words, the effect of spreading can 
be distributed between the couple of structure factors in an arbitrary manner. 

Likewise, the temporary energy f associated with [/("^(r) and contributing to 
(j54p can be presented in the form 

£(") = 1 /" p(r)t/(")(r)dr. (56) 
2 Jv 

According to (|5ip and ([5^ . one can see that formula ([55)) can also be rewritten as 

= 1 /" p(™)(r){/("-")(r)dr (57) 



at < TO < n. This fact justifies the complementary character of results ([5T|l and 
([52]). 

In terms of the space transformation it implies that two charged examples of 
transformed space interact either via (j57p or via the first term on the right-hand side 
of (|54p with account of ([55)1 . In the symmetric case of n = 2to both of these examples 
of space appear to be identical. From the standpoint of symmetry, such an event of 
the highest symmetry relative to the effect of spreading is the most beautiful. The 
corresponding symmetric case at n = 2 is the essence of the original treatment of 
Bertaut 

Nevertheless, the chief objective of spreading is to improve the calculation of 
electrostatic potentials in crystals. This is the reason that all the effect of spreading, 
notwithstanding is it single or multiple, should practically be attributed to the 
potential part of the Coulomb energy that is eventually in conjunction with the 
principal idea of Bertaut [27l [28] . 



5. General properties of bulk potentials at a multiple spreading 



The close connection between the Coulomb energy and electrostatic potentials results 
in the known fact that the potential at any point can be determined as a variational 
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derivative of the energy at hand with respect to the charge density at the same point 
[55] . With making use of relation ([55]) , it imphes that 

keeping in mind that this result is quite general and so it is numerically independent 
of the subscript n, as mentioned above. 

A similar result but with a distinct interpretation appears if we deal with the 
energy determined by formula ([57]) . In this case the revised version of ([55]) takes the 
form 

[/("-") (r), (59) 

where the restriction < m < n means that there are n + 1 different events 
associated with definition ([59]) . In other words, we have derived that the potential 
fields [7("~™)(r) with different superscripts can arise from a given 5^"). Note that 
along with the nth power of S{h) in ([39]) . this ambiguity for n = 2 was discussed 
earlier [52]. For completeness, it should be emphasized that the same potential field 
lj{n-m) ^^-^ can be also obtained from f '^"^ corresponding to different n. To this end, 
formula ([59)1 has to be rewritten as follows: 

[/("-'"^(r), (60) 



(5p("+'=)(r) 

where k > — m. Thus, issue ([59]) may be treated as a particular case of ([60| at fc = 0. 
Relations (l59)l and ([60]) enable one to render some debatable places associated with 
charge spreading more tractable. Indeed, according to ([M]) . each of the potentials 
occurring in ([55)1 and ([SO]) corresponds to the solution of Poisson's equation specified 
by the charge distribution /5("~™)(r) appropriate to the case. In this respect, these 
potentials are quite determinate. On the other hand, the connection between these 
potentials and the energies specified by ([57]) is also definite, despite the fact that 
different energies associated with the effect of charge spreading can be built up on the 
ground of the same potential field. This inference agrees with the conclusion known 
in the literature [22 [ [27l l28] . 

Interested in general spatial properties of potentials connected with the charge 
spreading, now we discuss the mean potential value defined as: 

U= - [ U{r)dr. (61) 

Substituting the first term on the right-hand side of p9|) into ([61]) . we encounter with 
the relation 

- / exp(27rift,r) dr = (SfiOj (62) 

V JccW 

where ShQ is the Kronecker delta. Formula (|62p is the fundamental relation of 
orthogonality describing the transformation from the real space representation to the 
reciprocal space one. Indeed, if ([3]) is substituted into ((4|), the result becomes the 
identity due to relation ([62)l . Hence, one can see that owing to the absence of the 
h = Q contribution to the first term on the right-hand side of ([55)) . this contribution 
does not affect the value of ([CTjl . 
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Considering the contribution of the second term on the right-hand side of ([35]) to 
(|6ip. we focus on the relation 



-f ^^Ef751-^*"^(l^^l)l (63) 



appearing in this case. Similar to the transformation from ([5]) to the integration 
over the unit cell along with the summation over i is transformed to the integration 
over all space again. As a result, expression ([55]) takes the form 

G(") = i / dr[i-/?(")(|r|)j, (64) 
V J Ir J 

where definition is taken into account and the corresponding shift r — > r — r' is 
suggested without changing the result. The integration over the angular variables of r 
in the first term in the square brackets is trivial, whereas in the second one it is readily 
performed if we go over to the new variabe defined by equation (IA.2P in [Appendix A[ 
Then we get 

G(") ^ 1^ ^ 



rdr - J <j{\ri\) . . . cr{\r„\) dri ...dr„ 



(i 



1 f'^ M 

- / r'^dr+ / rdr) , (65) 



where Q = |ri + . . . + r„|. If the first term in the square brackets in (|65p is formally 
multiplied by n integrals of form (j23p . then it can be combined with the second term 
therein. The result is as follows: 

G("' = ^ y a(|ri I)... a(|r„|)dri...dr„^'^(l-^) rdr. (66) 

The integration over r is straightforward here and we obtain 

= ^ /gVki|)...a(|r„|)dri...dr„. (67) 

It is significant that 

= + . . . + + 2(rir2) + . . . + 2(r„_ir„), (68) 

where all scalar products vanish after integrating over angles in (j67p , but all terms 
give equal contributions to (|67[) . Keeping relation (|23p in mind, we finally obtain 



3w 7 3u 
It is not surprising that result ([69]) looks like the mean potential of Bethe [29] addressed 
to n 'charge' distributions a{r) in a unit cell. It is important that G*^"^ turns out to 
be a constant. Thus after substituting ([63]) into formula ([39]) . result ([T]) arises and so 
this contribution is zero as well. The last term on the right-hand side of ([55]) is defined 
on a set of discrete points. Therefore its contribution to ([CT]) is of measure zero and 
so it is negligible. As a consequence, we deduce 

[7 = 0. (70) 

In other words, in uniform space zero mean charge, even with the effect of spreading, 
generates zero mean potential ^ i33j . 
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6. Simple exponential spreading 



For practical calculations some special representation of issue (|37|) is of interest. 
Indeed, as shown in [Appendix A[ formula ([37|) can be rewritten in the following 
recursion form 

n i<oo pR+r 

/2H(i?)= / a{r)rdr n^"-^\y)ydy, (71) 

where Q'^'^\R) — 1/R. Based on this relation, one can also obtain the limiting result 
useful in what follows: 

poo 

/2(")(o) = 47r / a{r)n^''-^\r)r'^dr. (72) 



In the particular case of n = 1 the values of W''^\R) and O'^^^O) follow from 
40| and ([TTjl and from ((72|) . respectively: 



H/(i)(i?) = 47r / a{r)r{r - R)dr, (73) 
Jr 

/•OO 

/2(i)(0) = 4tt a{r)rdr, (74) 



where equation (p3|) is employed. 

The values of W'-^^R) and /2(2)(0) are also of special interest. Their calculation 
is more tedious and is represented in [Appendix A[ The corresponding general results 
are as follows: 

/•oc />oo /^R pR—r 

W^^\R) = A-k^ / dr I dr'A{r,r')- dr dr'A{r,r') 



2 dr dr'B{r,r') , (75) 



iR+r 
/>oc />r 

n'^^\0) = 327r2/ a{r)rdr a{r'){r'f dr\ (76) 
Jo Jo 

where in formula (I75|) wc introduce the notations: 

j4(r, r') = a(r)a{r')rr' {R — r — r')^, (77) 

B(r, r') = a{r)a{r')rr'{R + r- r'f. (78) 

It is worth noting that J7'^^^(i?) determined by (l37|) and associated with M^(^)(i?) in 
form ([75|) through (j40|l . may be regarded as the energy of Coulomb interaction between 
two 'charge' distributions a{r) of the distance R apart [25]. Likewise, /2(^)(0) is 
appropriate to the energy of self-interaction, in agreement with the Bertaut treatment 

There is a large variety of spreading functions discussed in the literature 
[HliTlllSTllllllllliailTlllHllMllMlIMllST]. interested in principal aspects of charge 
spreading, here we first consider the simplest spreading function a{r) which falls off 
exponentially with the distance r — \r\ [571155] : 

Q,3 

'^(r) = -^exp{~ar), (79) 

providing that this function is normalized in compliance with (|25)) . Substituting ([75)1 
into pip , we readily derive 

5,.)^^l+(^)T^ (80) 
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As far as is concerned, we notice that this value is diniensionless in accord 

with its definition (I40p . It is then evident that this value can be cast in the form 

W^'^Hr) = W^"Hz), (81) 

where z is the diniensionless combination of R and the spreading parameter. In the 
present case it implies that z = aR. After inserting (|79|) into ((73| and ((74)l . we obtain 

W^W(z) = (l + |)exp(-z), (82) 

r2W(0) = |, (83) 

where result (|83|) follows from the combination of (|40|) and ([82|) as i? tends to zero. 
The case of ti = 1 arises upon substituting formulae (l80 |) -(f83 |) into (I39p . 

If the multiple spreading associated with rt = 2 is concerned, result ([80|) is still 
suitable. Substituting (1751) ™to expression ([751) . '^^ turn obtain 

(;,) = (i + ii£ + ^ + f!) exp(-z). (84) 



16 16 48 



The value of 



^^^^(i?) = ]^ - + + ^ + ^) -P(-^)- (85) 

is reconstructed from (|84l) with the help of (|40l) and ([84]) . It is important that formula 
(|85p describes the interaction between two identical exponential charge distributions of 
which centres are separated by the distance R, in accord with the inference mentioned 
above. In the limit of i? ^ formula ([85]) yields 



/2(2)(0) = ^, (86) 

in agreement with ([76l) . Based on relations ([80l) . ([84l) and ([86l) . expression ([39l) at 
n = 2 describes the potential of interest. 

Substituting ([55]) into formulae ([71]) and ([7^ . one can obtain the next generation 
of results appropriate to n = 3. They are of the form 

w ^ (' + H - w - w + IS + 4) '-(-^)' 

««'(0) = ||. ^ (88) 

where the transformation of n'^^^{R) to W''^'>{z) is carried out with making use of ([^O]) 
and ([8T|) again. Starting from equations ([71]) . ([72)1 and ([871) . results for n > 3 can be 
obtained in the same manner. 



7. Invariance of the Ewald approach 

Now we consider Gaussian functions of spreading. Interested in the spreading of the 
nth order, we introduce 

'^"W = (^)'^'exp(~Vr2) (89) 

that is normalized by condition ([23| . Substituting ([89]) into ([3T|) , integrating the result 
in Cartesian coordinates and keeping in mind the familiar Poisson integral [39j 
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we get 

S^{h)^exp( (91) 

With making use of (PT|) . relation takes the form 



that turns out to be independent of n. 

The consideration of n^'^^{R) in the particular cases of n = 1 and n = 2 may be 
performed basing on relations (|73p and (j75p . respectively. However, even at n = 2 
the corresponding relation is rather complicated. It is evident that the complexity 
will further enhance for n > 2. This obstacle can be overcome within the approach 
proposed by Boys [3D], where the integration over angular variables turns out to be 
much more efhcient if those variables are incorporated directly into the exponents of 
Gaussian functions, as shown in Appendix B 



Let us consider, for a moment, the spreading function without normalization: 

o'air) — exp(— ar^). (93) 



According to results (|B.4p and (|B.9p from Appendix B we then obtain 



2nl '■^ 



ni^\R)^ / exp{-ar^)dr, (94) 

a RJq 

^^\R)-^A-^%rM--y)<^r. (95) 



aV2\aJ RJq "^V 2 

Based on formulae ([55)1 and (jB.9p from [Appendix B| one can prove by induction 
that 

a^nl\ay J R Jq \ n / 

Now in (j96p we replace spreading function (|93p by the normalized one described by 
As a result, formula (j96p is transformed into 

n^''\R)^ r exp{-fi\^)dr. (97) 



Vt^R Jo 

It is clear that the dependence upon n disappears here. In the particular case of i? = 
relation ([57)1 yields 

/2(«)(o) = ^. (98) 

On the other hand, relation ([97|) can be identically rewritten as: 

^''(R) = ^-'^^, (99) 
where the complementary error function 

erfc(z) = ^/ exp(-w2Wy ^ noQ) 

just describes the value of W^("'(z) at z = /ii? due to the last equality that follows 
upon comparing (|99p with ([40[l and taking ([81]) into account. 
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If we substitute results 



and poop into then we derive 

2|»,|2 



f p{r') eTic{^i\R,\) dr' f 2^gj | 



(101) 



where Ji^ is defined by p6)) . On making use of formula (|10ip in (|53p . the specific 
energy takes the form 



27r^ |/i|2 



■ exp 



' /■ p(r)p(rO erfc(Aj|-R,|) 



drdr' -^Y.' 



(102) 



Equations (jlOip and (|102p are the classical formulae of Ewald |13J. We draw a 
conclusion that Gaussian spreading functions appear to be invariant with respect to 
their multiple application, without changing the functional form of the result. 

8. Multiple charge spreading in individual pair interactions 



According to ([38|) . the potential effect of spreading charges is separated from that of 
the initial ones in crystals. Therefore the final analytical results depend solely on a 
single value of n specifying the order of a multiple spreading at hand. The rate of 
convergence in dependence on n is a special subject that will be discussed elsewhere 

m- 

Here we concentrate our attention on another manifestation of charge spreading 
keeping in mind that in a single neutral object the charge spreading distribution 
may be treated as neutralizing a more compact charge of opposite sign |16j . Let us 
consider two complex objects of this sort, which are in general defined by the charge 
distributions 



pi(r)-Zi /('"i)(o,r) (o,r) 



(103) 



P2(r) 



(104) 



/(™-^)(o,r)-/("^)(o,r) 

where Zi and Z2 are the total charges describing either part of pi{r) and P2(''"), 
respectively, in accord with (|44p and (|46p . Here we restrict ourselves to a multiple 
application of a certain initial spreading function so that relations (|¥T|) - ([151) are 
taken into account in definitions (|103p and (|104p . The energy of interaction between 
charge densities (|103p and (|104p separated by the distance R can be written in the 
conventional form as 

Pi{ri)p2{r2) dridr2 



E{R) 



(105) 



l-R + ri -rsl 

On substituting (|103p and (|104p into (|105p and taking formula ([50]) into account, 
expression (jlOSp can be transformed into 

E{R) = ZiZ2W"'^+"'^\R) - /?(™i+"^)(i?) - l?("i+'"-^)(i?) 

Z1Z2 



/?("i+"^)(i?) 



R 



(106) 




Figure 1. The "interatomic" energy E, in units of aZiZ2, versus a dimensionless 
distance z = aR. Curve 1 is appropriate to equation JTOSj and exhibits the case 
of a simple exponential charge clouds neutralizing central point charges of the 
distance R apart. The case of the corresponding Gaussian charge clouds specified 
by formula II109I I is described by curve 2, providing that its minimum coincides 
with the minimum point of curve 1 at = 0.3586a that is marked by the vertical 
dotted line. The dashed lines indicate the points of intersection of these curves. 



Here each n^---\R) describes the interaction energy between the corresponding single 
terms in (|103[) and (|104[) and the transition to the last relation is performed by means 
of (j40|) . As a result, four different orders of spreading appear in this general case. Of 
course, a particular event at mi = m2 = and ni ~ n2 — I is of special interest. In 
this case W°{R) = due to (gOl), ^ and jSQl). Formula (fTOe]) is then converted into 

Z1Z2 



E{R) 



R 



2W^^\R)-W^^\R) 



(107) 



where only two consecutive orders of spreading happen. 

The particular case associated with a simple exponential spreading arises after 
inserting relations (|82|) and (|84p into (|107|) . Keeping (|8T|) in mind, we then obtain 



exp(-z). 



(108) 



5z 3z2 
^ 16 ^ Te" ~ 48J 

where z = aR. Likewise, the case appropriate to Gaussian spreading functions arises 
upon substituting (jlOOp into (|107p . with taking underlying definition ([55]) into account. 
Starting from o-i(r) in ((89)) . we readily obtain the following result 



Eg{z) - 



2erfc(z) -erfc(z/V2) 



(109) 



where z = /ii? and 2(/i/V2)^ in definition ([55)1 so as to describe the last term 

in the parentheses in (|109p . 

Comparing issues (|108p and (|109p , we recognize one more universality, which now 
corresponds to the energy described by the simple exponential spreading, where the 
same exponent turns out to be typical of both the terms addressed to W^^^R) and 
W^"^^ (R). Actually, it is not surprising because the next sample of this set, i.e. W'^^^ (i?) 
described by ([57]) . is specified by the same exponential decrease. Conversely, if the 
energy is determined by Gaussian spreading functions, then either of the contributions 
to (|109p is specified by its own predominant law of decrease. 




Figure 2. The parameter K defined by formula (IllOI I versus a dimensionless 
distance z = aR. The position of energy minima, where K has a local minimum, 
is pointed out by the vertical dotted line. 



It seems to be interesting to recall that there is sometimes a tendency towards 
changing a simple exponential electron spreading by a Gaussian one in molecular 
calculations, where the contribution of the exchange interaction can then be evaluated 
in a much simpler manner [40l HH [42j [43] . Although the direct Coulomb interaction is 
still the subject of our interest, this is the reason to compare results (|108|) and (|109p in 
more amount of detail. To this end, we plot the corresponding energy curves together, 
as shown in figure [TJ The shape of either of these curves is quite natural. Indeed, the 
energy is positive and its value tends to infinity as R drops to zero. On the other hand, 
if R grows, then the energy eventually becomes negative, attains at its minimum value 
and farther falls off to zero in magnitude, being still negative. The latter is a direct 
consequence of the fact that just in general relation (|107p the contribution of W^^\R) 
as a function of R is always more diffuse than that of W^^\R). In order to compare 
both the curves, we have shifted a minimum point of the curve describing quantity 
(|109P to the value of i? = 3.745/a. We see that a minimum of curve 2 corresponding 
to Eq™ = —0.01715aZiZ2 is much deeper than that of curve 1 with the value of 
E^'"" = -0.0098050:^1^2 and this effect is described by the ratio 1.749. Moreover, 
the curvature of curve 2 at its minimum point is also greater than that of curve 1. As 
a result, there are two points of intersection between those curves which take place at 
R = 2.781/a and at i? = 6.830/q;, with the energy values E = -0.000652aZiZ2 and 
E = — 0.00194q;ZiZ2, respectively, as shown in figured] as well. 

Of course, the energy Eq{z) can be further scaled by the factor ^'"""/i?™'" so 
as to simulate the behaviour of Es{z). The comparison of both these energies is then 
specified by the relative value of the form 

W |_ t^V /^min sv ^J|£;min| |£;mm| |£;mm| ^ ^ 

The behaviour of K{z) as a function of z is shown in figure|2] We see that K{z) turns 
out not to be monotonic in the vicinity of z corresponding to minima of the energies at 
hand. This fact can be important upon minimizing a total energy modified by other 
energy contributions. In this case the effect driven by simple exponential spreading 
functions and that driven by Gaussian spreading functions are expected to be far from 
being proportional. 
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9. Conclusion 



In summary, it is shown that the effect of charge spreading proposed by Bertaut 
[2T| [27] can be utihzed in a multiple manner. It means that the problem how many 
times a given spreading function is applied to the original charge distribution in a 
crystal is not of principle. Nevertheless, the tendency towards increasing the rate of 
convergence upon multiple charge spreading is just recognized by Bertaut [27] and will 
be confirmed elsewhere [32]. This result is not trivial. Presumably, it is associated 
with an idea that there is an optimum spreading configuration with very diffuse tails. 
In this connection, the fact that all the effects driven by a Gaussian spreading function 
are reproduced in the same functional form, regardless of its multiple application, may 
anyhow point to an optimum character of a Gaussian spreading. 

Here we also recognize that a certain spreading, either single or multiple, may be 
attributed to every point of space and so may be treated as a uniform transformation 
of space. It is significant that the general relation between electrostatic potentials and 
specific Coulomb energies in crystals, as well as zero value of the mean potential there, 
turns out to be invariant with respect to such a transformation. 

It is evident that the application of a multiple charge spreading to the problem of 
lattice summation is nothing but a fruitful approach to that problem. It implies that 
the final results of lattice summation are to be independent of the shape of spreading. 
However, it is not the case if the charge spreading is regarded as a real property of at 
least a pair of complex neutral physical objects connected by the Coulomb interaction. 
In this event the replacement of a natural shape of, for example, an electron cloud 
with a more artificial shape must be performed with caution. 



Appendix A. Some general relations for charge spreading 

Equation ([37]) can be readily rewritten as 

n^''\\R\)^ j a(|r|)/?("-i)(|H + r|)dr. (A.l) 
With making use of spherical coordinates of r, we assume that 

\R + r\^ {R^ + r^ + 2Rr coney = y (A.2) 

and go over from the variable to a new variable y. The result of integration over y 
is then of form fTl]) . On the other hand, if il = 0, then the integration over angular 
variables of r in (|A.ip is trivial and we obtain issue ([7^ . 

The case of fl'^^^ (R) is straightforward and is described by 

n^^\R) ^ ^ - Att J <T(r)(-^ - l)rdr. (A.3) 

Equations ([75]) and ([73]) follow therefrom. If we are interested in J7'^^^(i?), then the 
employment of (|A.3[) in the general relation ([7T]) gives rise to 

Q^^\R)^ Q^^\R)-^ I a{r)J{R,r)rdr, (A.4) 



R 

where we utilized definition (jA.3P again and 

pR+r pOQ 

J{R,r)= dy a{r'){r' -y)r' dr'. (A.5) 

J\R-r\ Jy 
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If we interchange the order of integration over r' and y here, then the integration over 
y is straightforward and we obtain 

j(R^r) = - a{r'){\R^r\-r'fr'dr'~- a{r'){R + r - r'fr' dr' . (A.6) 

2 J\R-r\ 2 jB,+r 

In turn, based on (|23p . it is expedient to rewrite expression (jA.3|) in an identical form 

1 1 (^™-2 /'OO /•oo 

n(^){R) = — a{r)rdr a{r'y [r{r' - R)\dr' . (A.7) 

^ Ip JR 

Substituting (jA.6[) and (jA.7|l into (IA.4P and combining the integral terms, we arrive 
at the result in the most symmetric form given by formulae ([75)1 . ([77)1 and ([75]) . in 
agreement with (|40l) . On the other hand, relation ([75]) for J?'^^^(0) appears directly 
upon substituting (jA.3p into ([7^ . 

Appendix B. Coulomb interaction between Gaussian functions 

Here we follow the treatment of Boys IJD] • Let us consider two charge distributions 

cra(r) = exp(-Q!r^), cr/3(r) = exp(-/3r^). (B.l) 
According to ([57)) . they determine the values 

= j cj^{r)nf \\R + r\)dr, (B.3) 

where i? = and r — \r\. 

The particular form of the distributions in (jB.ip enables one to go over to the 
variable r' = r + i? in equation (|B.2p . The integration over the angular coordinates 
of r' is then trivial there and wc obtain 

27r 

^fP^^)^J^j^ exp(-/?x2)dx, (B.4) 

where x = \r'\. Inserting (jB.4p into (|B.3p and operating further in the same manner, 
we derive 



where 



/oo ry 
exp[-~a{R~yf]dy / exp(/3x^) dx. (B.6) 
oo Jo 

l|B.6p with respect to R, then we obtain 

f°° d ( ~i /"^ 

j —[exp[-a{R-yf]\dy exp{Px^)dx. (B.7) 



If we differentiate equation (|B.6p with respect to R, then we obtain 
dI{R) r°° , ^ . rv 

dR ' 

Integrating the right-hand side of (jB.7p by parts, we easily reach 



dI(R) I TT / a/3 

— ' ; exp ' 



dR y a + f3 \ a + P 
Note that 7(0) = follows from (|B.6p and specifies the further integration of (|B 



with respect to R. On inserting the result of integration into equation (jB.Sp . the final 
issue takes the form 



n^)^R) ^ / expf-^) dz (B.9) 



a/Sy^a + l3R Jo ^ a + /3/ 
that is naturally symmetric with respect to a and /3. 
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